Analytically determining frequency and amplitude of spontaneous alpha oscillation in Jansenʼs neural mass model using the describing function method
Xu Yao1, 2, Zhang Chun-Hui1, Niebur Ernst3, Wang Jun-Song1, †
School of Biomedical Engineering, Tianjin Medical University, Tianjin 300070, China
Qingdao Stomatological Hospital, Department of Medical Technology Equipment, Qingdao 266001, China
Zanvyl Krieger Mind/Brain Institute and Solomon Snyder Department of Neuroscience, Johns Hopkins University, Baltimore 21218, MD, USA

 

† Corresponding author. E-mail: wjsong2004@126.com

Abstract

Spontaneous alpha oscillations are a ubiquitous phenomenon in the brain and play a key role in neural information processing and various cognitive functions. Jansenʼs neural mass model (NMM) was initially proposed to study the origin of alpha oscillations. Most of previous studies of the spontaneous alpha oscillations in the NMM were conducted using numerical methods. In this study, we aim to propose an analytical approach using the describing function method to elucidate the spontaneous alpha oscillation mechanism in the NMM. First, the sigmoid nonlinear function in the NMM is approximated by its describing function, allowing us to reformulate the NMM and derive its standard form composed of one nonlinear part and one linear part. Second, by conducting a theoretical analysis, we can assess whether or not the spontaneous alpha oscillation would occur in the NMM and, furthermore, accurately determine its amplitude and frequency. The results reveal analytically that the interaction between linearity and nonlinearity of the NMM plays a key role in generating the spontaneous alpha oscillations. Furthermore, strong nonlinearity and large linear strength are required to generate the spontaneous alpha oscillations.

1. Introduction

Alpha oscillation, defined as periodic fluctuations within the alpha band (8 Hz–12 Hz),[1] is ubiquitous in the brains of animals and humans.[2,3] Spontaneous alpha oscillations are observed in waking subjects and are most prominent in subjects that are relaxed with their eyes closed and otherwise mentally inactive. Since the first human electroencephalogram (EEG) recordings were performed by Hans Berger more than 90 years ago,[1] many studies have tried to establish correlations between the characteristics of alpha oscillations and cognitive processes.[419] Although several studies have shown that alpha oscillations are related to a wide range of basic and higher cognitive processes and play a functional role in both behavioral and neuronal responses, the mechanism of alpha oscillation generation is not fully understood.

Computational studies[2025] are effective for elucidating the mechanism by which neural oscillations occur due to their ability to model complex neurological phenomena with relative ease.[26] In this context, many biophysical models have been developed to understand the generation of neural oscillation.[2732] As a typical computational neural model, the neural mass model (NMM) was initially proposed to study the origin of the alpha rhythm based on a biologically plausible parameterization of the dynamic behavior of the layered neocortex.[33,34] Subsequently, the NMM has been improved and extended to model various cortical electrical activities.[35] Jansenʼs NMM is an extended version of the neural mass model that adds an excitatory feedback loop in parallel to the inhibitory feedback loop in the original version of the NMM.[36]

To date, Jansenʼs NMM (for simplicity of denotation, we will use NMM to refer to Jansenʼs NMM for the remainder of this paper) has been extensively used to study alpha oscillation mainly through numerical study, such as resonance and bifurcation analyses.[27] Resonance analyses are conducted by linearizing the sigmoid function in the NMM.[35] However, this produces a stimulus-induced alpha response and not a spontaneous oscillation, which is dominated by linear dynamics and reflect fluctuations about an equilibrium-state.[37] Most of bifurcation studies use the input of the NMM as a bifurcation parameter to study neural oscillation; therefore, the responses produced in these studies are also stimulus-dependent oscillations. Spontaneous alpha oscillation is related to limit cycle dynamics[27,28] and is primarily determined by the intrinsic properties of the NMM. Nonlinearity is the fundamental reason for the occurrence of limit cycles, however, it is not fully understood how the nonlinear part of the NMM interacts with its linear part to cause spontaneous alpha oscillations.

In the present study, we use the describing function method to analyse the spontaneous alpha oscillations in the NMM. Our main contribution is to propose an analytical method to accurately determine frequency and amplitude of the spontaneous neural oscillation. To do so, we take Jansen–Ritʼs neural mass model as an example to demonstrate the method. Our results suggest that the interaction between linearity and nonlinearity of the NMM plays a key role in generating spontaneous alpha oscillations and that, furthermore, strong nonlinearities and large linear strength are required. These results provide an insight into the mechanism of spontaneous alpha oscillations and the role of the linear and nonlinear parts in generating them. Additionally, these results provide testable hypotheses for future experimental studies.

2. Methods
2.1. Model

Jansenʼs neural mass model describes a small local network representing a cortical region. As illustrated in Fig. 1,[36] the model consists of three interacting neural populations: a feed-forward loop (the middle part), an excitatory feedback loop (the top part), and an inhibitory feedback loop (the bottom part). Jansen–Ritʼs neural mass model (NMM) shows that the activity of the excitatory population in the feedforward loop is regulated by the interaction between the interlinked excitatory and inhibitory feedback loops. The connectivity constants C1, C2, C3, and C4 denote the interactions between the feed-forward and feedback populations in the NMM and account for the average number of synaptic contacts. The input p(t) is an excitatory input from neighboring populations and is modelled as Gaussian white noise. The output y(t) represents the average synaptic activity of the pyramidal cells, which is explained as an EEG signal.[36]

Fig. 1. Schematic diagram of the neural mass model.

In the NMM, each population receives input in the form of an average postsynaptic membrane potential from the other neural populations and converts this membrane potential into an average firing rate. Therefore, each population is composed of the following two parts: an excitatory or inhibitory linear dynamic function or , respectively] and a nonlinear static sigmoid function S(v). The former transforms the pre-synaptic mean firing rate into the post-synaptic mean membrane potential and the latter transforms the mean membrane potential into the mean firing rate.

The impulse response functions of the excitatory and the inhibitory dynamic functions are given by the following equation[36] where u(t) is the Heaviside function, He and Hi denote the excitatory and inhibitory average synaptic gains, respectively, and τe and τi are the lumped representations of the sums of the membrane average time constants.

The sigmoid function S(v) is described by the following equation:[38] where e0 is the maximum firing rate, v0 is the post-synaptic potential corresponding to a firing rate of e0, r is the steepness of the S(v) function, and v is the average pre-synaptic membrane potential.

The parameter values of the NMM can be estimated from the empirical data,[34,36,39,40] and are shown in Table 1.[36] The focus of this study is to investigate spontaneous oscillations and not stimulus-dependent oscillations; therefore, the input p(t) is set to be zero. For further details of the NMM see Ref. [36].

Table 1.

Physiological explanation and standard values of NMM parameters.

.
2.2. Describing function method

The describing function method is a nonlinear analysis method that is widely used to analyse spontaneous oscillations of nonlinear dynamical systems.[41] Recently, this method was used to analyse the mechanism of action of deep brain stimulation,[42] the behavior of a seizure neural mass model[43] and a model of deep brain stimulation.[44] The main principle of the describing function method is that in response to a sinusoidal input, most nonlinear system will produce a periodic signal (which is not necessarily sinusoidal) with the frequency being the harmonics of the input frequency. Therefore, the describing function can be viewed as an extension of the frequency response to nonlinearity, which can be approximated by a frequency-dependent equivalent gain.

Supposing that the input to a nonlinear function is , the output q(t) will be periodic with a fundamental period equal to that of the input. Consequently, the Fourier series is described as follows: where is the direct current component, is the n-th harmonic component of q(t), and and are the coefficients of the Fourier transformation.

Because the linear part of a dynamical system is generally of low-pass, the high frequency component of q(t) decreases rapidly as the frequency increases. Thus the output q(t) of the nonlinear element can be approximated by the first fundamental component of this series and calculated as follows:[45] The describing function is defined as the ratio of the amplitude of the first fundamental component of the output of the nonlinear element to the amplitude of the sinusoidal input signal.[45] Its mathematical formulation is expressed as follows: More details of the describing function method can be found in Ref. [45].

Within the theoretical framework of the describing function method, the block diagram of a nonlinear system is shown in Fig. 2, where N(A) is the describing function of the nonlinear part in the nonlinear system, and L(s) is the transfer function of the linear part in the nonlinear system.

Fig. 2. (color online) Block diagram of a nonlinear system.

The characteristic equation of a nonlinear system (shown as in Fig. 2) is as follows: Letting and substituting it into Eq. (6), the characteristic equation of a nonlinear system in the frequency domain is obtained as follows: According to the Nyquist criterion,[41] the solution of Eq. (7) corresponds to the limit of the cycle oscillations; therefore, solving this equation, the amplitude and frequency of the spontaneous oscillation are obtained.

2.3. Simulation protocol

All simulation analyses were conducted using MATLAB software (www.mathworks.com). Acording to Fig. 1, we used the Simulink toolbox of MATLAB to construct the simulation model of the NMM, in the form of a block diagram. In the simulated model, the excitatory and inhibitory synapse models are implemented by Eq. (8), and the Sigmoid functions by Eq. (2).

The simulation parameters of the NMM were set and described in Table 1. To analyze and simulate the spontaneous alpha oscillations in the NMM, the input of the NMM was set to be zero. To trigger the spontaneous oscillation, in the initial state the input was an impulse with arbitrarily small amplitude and width. The simulated model was solved using fourth order Runge–Kutta method with fifth order error control (Matlab function ode45).

We used a variable-size simulation time-step with maxim step size 0.001 s. Initial conditions were set to be zero.

Nyquist diagrams were plotted using the MATLAB function “nyquist”.

3. Results

In this section, we use the describing function method to investigate the mechanisms underlying spontaneous alpha oscillations in the NMM. First, we reformulate the NMM (shown in Fig. 1) by approximating the sigmoid function using its describing function. Second, we conduct a theoretical analysis to determine whether spontaneous alpha oscillations occur and, furthermore, to determine the amplitude and frequency of the oscillation.

3.1. Frequency domain description of the NMM using describing function method
3.1.1. Transfer function of the linear part in the NMM

By conducting the Laplace transform with Eq. (1), we derive the transfer function, s domain model, of the excitatory and inhibitory synapses as follows:

3.1.2. Describing function of the sigmoid function in NMM

The sigmoid function S(v) in the NMM is odd symmetric; therefore, the coefficients of all the even parts are equal to zero, which means that and equation (5) can be simplified into Furthermore, we can obtain the describing function of S(v) as follows: The sigmoid function S(v) in the NMM can be written using the following form: Substituting Eq. (11) into Eq. (10) yields the following equation: According to Eq. (12), we can draw the curve of the describing function N(A) of the sigmoid function S(v) as shown in Fig. 3.

Fig. 3. (color online) Describing function N(A) of the sigmoid function S(v)).
3.1.3. Reformulating NMM in frequency domain using describing function method

To use the describing function method to analyze spontaneous alpha oscillations in the NMM, we need to transform the NMM into the standard form of a nonlinear system (as shown in Fig. 2).

First, we use N(A) to replace the three sigmoid functions S(v) in the NMM, which transforms Fig. 1 into Fig. 4(a).

Fig. 4. (color online) Block diagrams of the reformulated NMM using the describing function method. Panel (a) is a reformulated form of the NMM shown in Fig. 1 and is derived by replacing the nonlinear function S(v) using its describing function N(A). Panel (a) is further transformed into panel (b) by approximating the two describing functions N(A) in the excitatory and the inhibitory feedback loops. Panel (c) is an equivalently transformed version of panel (b), where G(s) is the transfer function of the linear part of the NMM (blue dashed rectangle) and N(A) is the nonlinear part of the NMM.

Second, we reformulate the NMM (shown in Fig. 4(a)) into the standard form of a nonlinear system (illustrated in Fig. 22), which is composed of one linear unit and one nonlinear unit. Note that the output y(t) of the NMM is directly linked to the sigmoid function S(v) in the feed-forward loop of the NMM. Therefore, the sigmoid function S(v) in the feed-forward loop is the nonlinear unit and the other functions in Fig. 4(a) correspond to the equivalent linear unit. We use bifurcation analysis to determine the amplitude of signal in Fig. 4(a) and obtain . In Fig. 4(a), the inputs of N(A) in the excitatory and the inhibitory feedback loops are and , respectively; therefore, their amplitude values are and , respectively. According to Eq. (12), we approximate the describing function N(A) in the excitatory and the inhibitory feedback loops in Fig. 4(a) as and , respectively. Therefore, we derive the approximated form of Fig. 4(a) (shown in Fig. 4(b)), which is composed of the nonlinear unit N(A) (red part) and the linear unit (blue part).

Third, we equivalently transform Fig. 4(b) into Fig. 4(c), which is the block diagram of the approximated NMM using the describing function. In Fig. 4(c), G(s) is the transfer function of the equivalent linear unit in the NMM and N(A) is the describing function of the nonlinear part in the NMM. According to Fig. 4(c), we can derive G(s) as follows:

3.2. Analytical results

In this section, the Nyquist criterion can be extended to understand the behavior of the NMM whose nonlinearities have been approximated using the describing function.[45] According to Fig. 4(c), we derive the characteristic equation of the reformulated NMM using the describing function method as follows: Letting and substituting it into Eq. (14), the characteristic equation of the NMM in the frequency domain is obtained as follows: where According to the Nyquist criterion,[41] the solution of Eq. (16) corresponds to the limit of the cycle oscillation (i.e., the spontaneous oscillation), and it yields the amplitude and the frequency of the spontaneous oscillation.

Usually, it is difficult to directly solve Eq. (16). We can draw the curve and in the complex plane according to Eqs. (12) and (13), respectively. Note that is real; therefore, always lies on the real axis. The intersection point between the curve and the line , in the form of , corresponds to the spontaneous oscillation in the NMM. Therefore, we can further derive The analytical results are shown in Fig. 5(a). The blue curve is the Nyquist diagram of the linear part , and the red line is the negative reciprocal of the describing function , which is a straight line that is coincident with the negative real axis and is parameterized as a function of amplitude A of the input signal. It is observed that the curve and the line intersect at the point . Substituting into Eq. (17) yields Therefore, solving Eq. (18), we can determine that the frequency and the amplitude of the spontaneous oscillation in the NMM are /s (i.e., and , respectively.

Fig. 5. (color online) Analytical results of the spontaneous alpha oscillations in the NMM using the describing function method. (a) Analytical results. The blue curve denotes the Nyquist diagram of the linear part , where the array represents the increase of angular frequency ω. The red line refers to the negative reciprocal of the describing function of sigmoid function S(v), which is calculated as . Intersection between the blue curve and the red line corresponds to the spontaneous oscillation in the NMM. (b) Simulation results of the spontaneous alpha oscillations in the NMM. (c) A comparison between analytical and simulation results.
3.3. Simulation results

To verify the validity of the analytical results, we conduct time–domain simulations of the NMM using MATLAB software (for details see Section 2). The simulation results are shown in Fig. 5(b), which illustrates that the frequency and amplitude of the spontaneous oscillation are f = 8.62 Hz and A = 10.41 mV, respectively.

A comparison between the analytical results and the simulation results is shown in Fig. 5(c), which demonstrates that the analytical results are in agreement with the simulation results.

3.4. Effect of linear strength and degree of nonlinearity

Figure 5(a) shows that the interaction between the linear part and nonlinear part of the NMM plays a critical role in generating the spontaneous alpha oscillation. According to Fig. 3 and Eq. (2), with the decreasing of degree of nonlinearity of the sigmoid function, the value of will become smaller. As shown in Fig. 6, this means that the red line will move to the left, which means that there is no intersection between the red line and the blue curve. This reveals that a strong nonlinearity is required to generate spontaneous alpha oscillations. According to Fig. 5(a) and Eq. (1), with the decreasing of strength of the excitatory or the inhibitory linear dynamic functions in Eq. (1), the value of the amplitude of will become smaller. As shown in Fig. 7, this means that the blue curve will contract, with the result that there is no intersection between the red line and the blue curve. This shows that large linear strength is also required to generate spontaneous alpha oscillations.

Fig. 6. (color online) Effect of the degree of nonlinearity. We set r = 0.2, smaller than its normal value r = 0.56, representing a decreased nonlinearity of the sigmoid function. In this case, there is no intersection between the red line and the blue curve. Inset: the output of the NMM.
Fig. 7. (color online) Effect of the linear strength. We set , smaller than its normal value , representing the decrease of the linear strength. In this case, there is no intersection between the red line and the blue curve. Inset: the output of the NMM.
3.5. Effect of model approximating error

In this study, we use the describing function N(A) to approximate the sigmoid function S(v) in the excitatory and the inhibitory feedback loops. We estimate the amplitude of in Fig. 1 and further determine the amplituded of the input signal of the sigmoid function S(v) in the two feedback loops. Therefore, the estimated amplitude value of will exert an influence on the approximated N(A) in the two feedback loops. Furthermore, the estimated amplitude value of has an effect on the results of the oscillation analysis of the NMM.

In this section, we discuss the effect of the model approximation error on the spontaneous oscillations in the NMM. To quantitatively evaluate the analytical precision, we define the analytical precision measure as , where and f are the analytical and actual values of the oscillation frequency, respectively.

The analysis results are shown in Fig. 8, where Figure 8(c) shows the effect of the scale of variation in the estimated amplitudes of on the analytical results. These results show that the analysis results are accurate despite some model approximation errors; therefore, the analysis scheme employed in this study is robust against the estimated errors of .

Fig. 8. (color online) Effects of the model approximation errors on spontaneous oscillation in the NMM. (a) Analysis of the spontaneous oscillations in NMM with different model approximation errors of y0. (b) Frequencies of the spontaneous oscillations with different model approximation errors of y0. (c) Relative errors of the frequencies of the spontaneous oscillations with different model approximation errors of y0, where , and is the estimated amplitude of y0.

Figure 8(c) further shows that the analytical result of the oscillation frequency in the case that the estimated amplitude of is larger than the actual value, is more accurate than the result in the case that the estimated amplitudes of is smaller than the actual value. This is, as shown in Fig. 3, because that the describing function of the Sigmoid function with large amplitude input is more robust against the variation in the estimated amplitudes of than that with small amplitude input.

4. Discussion

In the present study, we take a neural mass model as a test bed to demonstrate how to apply the describing function method to the analysis of the spontaneous alpha oscillations and determine their frequency and amplitude. The proposed method can be easily extended to other neuronal circuits.

For a large scale feedback neuronal circuit,[46] there often exists a delay in the signal pathway due to dynamical relaying[47] and signal propagation.[48,49] We can incorporate the delay into the transfer function of the linear part as shown in Fig. 2, and employ the describing function method to elucidate the role which it plays in generating the neural oscillations and regulating the amplitude and frequency. We predict that the delay enhances the generation of neural oscillations, since it reduces the linear stability of the feedback system, and that it plays a crucial role in quantitatively regulating the amplitude and frequency, since a delay exerts an important influence on the Nyquist plot of the open loop transfer function. We will present these results in a future study.

Autaptic connection is a typical structural pattern, and provides a self-feedback loop for neuronal circuits.[50,51] Previous studies have revealed that autaptic connections have an important influence on neuronal dynamics[50,51] and information processing.[52] Within the describing function theoretical framework, we can analytically elucidate the role that autaptic connections play by combining it into the linear part of the feedback system.

As many studies demonstrated, normal oscillations provide a substrate for neuronal information processing, such as information transmission[53] and neuronal computation.[54] On the other hand, abnormal oscillations are related to pathological diseases, such as epilepsy[3,25,5559] and Parkinsonʼs disease.[60] As is well known, the plasticity in neuronal circuits plays a central role in maintaining brain functions. Due to synaptic plasticity,[61] structural plasticity, and intrinsic plasticity, respectively, the parameters in the feedback loops of neuronal circuits, such as the coupling strength and connectivity parameters in the linear part, and the parameters in the sigmoid nonlinear unit, are often time-varying. The present study reveals that the describing function method can build a relationship between neural oscillations and the parameters in the feedback pathways of neuronal circuits including excitatory and inhibitory feedback, and self-feedback in the form of autaptic connections.[50,51] Thus the proposed method may bridge between these neural plasticity mechanisms of brain function and neuronal disease. In a disease state, the neural system may lose its stability and pathological oscillation results.[54,59,62] The proposed method can determine the transition condition from the stable state to the unstable state and the oscillation characteristics. Thus, the method has the potential to understand the dynamical transformation of brain neural activities, such as from the normal state to disease states.

5. Results

In this study, we use the describing function method to analyse spontaneous alpha oscillations in the NMM. The results demonstrate that the describing function method is useful in determining whether or not spontaneous alpha oscillations occur and in accurately determining the amplitude and frequency of the oscillations. The results also highlight that the interaction between the linear part and nonlinear part of the NMM plays a critical role in generating the spontaneous alpha oscillations. Furthermore, strong nonlinearity and large linear strength are required to generate the spontaneous alpha oscillations.

This proposed analytical method can build the relationship between the frequency and amplitude of the spontaneous oscillation and the model parameters in a more transparent manner than the numerical methods employed in most computational studies. Therefore this approach can provide an insight into how the model parameters exert effects on the frequency and amplitude of the spontaneous neural oscillation. It thus has the potential to gain an in-depth insight into the mechanism under the spontaneous oscillations. To sum up, this analytical method provides a novel theoretical framework which can be extended to the study of spontaneous oscillations in other computational neural models. We will develop more examples in future work.

Reference
[1] Berger H 1929 European Archives of Psychiatry and Clinical Neuroscience 87 527
[2] Klimesch W 2012 Trends Cogn. Sci. 16 606
[3] Wang J S Wang M L Li X L Niebur E 2015 Chin. Phys. 24 038701
[4] Jensen O Bonnefond M 2013 Trends Cogn. Sci. 17 10
[5] Samaha J Bauer P Cimaroli S Postle B R 2015 Proc. Natl. Acad. Sci. 112 8439
[6] Haegens S Händel B F Jensen O 2011 J. Neurosci. 31 5197
[7] Sauseng P 2012 Current Biology 22 R306
[8] Payne L Guillory S Sekuler R 2013 Journal of Cognitive Neuroscience 25 1463
[9] Brinkman L Stolk A Dijkerman H C de Lange F P Toni I 2014 J. Neurosci. 34 14783
[10] Sadaghiani S Scheeringa R Lehongre K Morillon B Giraud A L D’Esposito M Kleinschmidt A 2012 J. Neurosci. 32 14305
[11] Haegens S Nácher V Luna R Romo R Jensen O 2011 Proc. Natl. Acad. Sci. 108 19377
[12] Neuling T Rach S Wagner S Wolters C H Herrmann C S 2012 Neuroimage 63 771
[13] Palva S Palva J M 2007 Trends in Neurosciences 30 150
[14] Foxe J J Snyder A C 2011 Frontiers in Psychology 2 154
[15] Buchholz V N Jensen O Medendorp W P 2014 Frontiers in Human Neuroscience 8 446
[16] Lou B Li Y Philiastides M G Sajda P 2014 NeuroImage 87 242
[17] Strauss A Kotz S A Scharinger M Obleser J 2014 Neuroimage 97 387
[18] Spaak E deLange F P Jensen O 2014 J. Neurosci. 34 3536
[19] Roux F Uhlhaas P J 2014 Trends Cogn. Sci. 18 16
[20] Hindriks R vanPutten M J 2013 Neuroimage 70 150
[21] Hindriks R Woolrich M Luckhoo H Joensson M Mohseni H Kringelbach M Deco G 2015 NeuroImage 106 328
[22] Becker R Knock S Ritter P Jirsa V 2015 PLoS Comput. Biol. 11 e1004352
[23] Mayhew S D Ostwald D Porcaro C Bagshaw A P 2013 Neuroimage 76 362
[24] Bhattacharya B S Coyle D Maguire L P 2011 Neural Netw. 24 631
[25] Wang J Niebur E Hu J Li X 2016 Sci. Rep. 6 27344
[26] Gerstner W Sprekeler H Deco G 2012 Science 338 60
[27] Grimbert F Faugeras O 2006 Neural Comput. 18 3052
[28] Spiegler A Kiebel S J Atay F M Knösche T R 2010 Neuroimage 52 1041
[29] Goodfellow M Schindler K Baier G 2012 NeuroImage 59 2644
[30] Nevado-Holgado A J Marten F Richardson M P Terry J R 2012 Neuroimage 59 2374
[31] Touboul J Wendling F Chauvel P Faugeras O 2011 Neural Comput. 23 3232
[32] Bhattacharya B S Cakir Y Serap-Sengor N Maguire L Coyle D 2013 Neurocomputing 115 11
[33] Da Silva F L Van Lierop T Schrijer C Van Leeuwen W S 1973 Electroencephalogr. Clin. Neurophysiol. 35 627
[34] Da Silva F L Hoeks A Smits H Zetterberg L 1974 Kybernetik 15 27
[35] Babajani-Feremi A Soltanian-Zadeh H 2010 Neuroimage 52 793
[36] Jansen B H Rit V G 1995 Biol. Cybern. 73 357
[37] Hindriks R Bijma F van Dijk B van der Werf Y D van Someren E J van der Vaart A 2011 Neuroimage 57 440
[38] Ursino M Cona F Zavaglia M 2010 NeuroImage 52 1080
[39] Jansen B H Zouridakis G Brandt M E 1993 Biol. Cybern. 68 275
[40] van Rotterdam A Da Silva F L Van den Ende J Vier gever M Hermans A 1982 Bull. Math. Biol. 44 283
[41] Franklin G F Powell J D Emami-Naeini A 1994 Feedback control of dynamics systems 6 New Jersey Pearson Higher Education 157 175
[42] DePaor A M Lowery M M 2009 Ieee. T. Bio-med. Eng. 56 2717
[43] Shayegh F Bellanger J J Sadri S Amirfattahi R Ansari-Asl K Senhadji L 2013 Journal of Medical Signals and Sensors 3 2
[44] Davidson C M dePaor A M Lowery M M 2014 IEEE T. Bio-med. Eng. 61 957
[45] Hu S S 2007 Automatic Control Principle 5 Beijing Science Press in Chinese
[46] Siegel M Donner T H Engel A K 2012 Nat. Rev. Neurosci. 13 121
[47] Vicente R Gollo L L Mirasso C R Fischer I Pipa G 2008 Proc. Natl. Acad. Sci. USA 105 17157
[48] Guo D Q Wang Q Y Perc M 2012 Phys. Rev. 85 061905
[49] Chacron M J Longtin A Maler L 2005 Phys. Rev. 72 051917
[50] Guo D Q Wu S D Chen M M Matjaž P Zhang Y S Ma J L Cui Y Xu P Xia Y Yao D Z 2016 Sci. Rep. 6 26096
[51] Guo D Q Chen M M Perc M Wu S D Xia C Zhang Y S Xu P Xia Y Yao D Z 2016 Europhys. Lett. 114 30001
[52] Guo X Yu H Wang J Liu J Cao Y Deng B 2017 Physica A Statistical Mechanics & Its Applications 482 308
[53] Buzsáki G Schomburg E W 2015 Nat. Neurosci. 18 484
[54] TP V R K LF A 2005 Ann. Rev. Neurosci. 28 357
[55] Chen M M Guo D Li M Ma T Wu S Ma J Cui Y Xia Y Xu P Yao D 2015 PLoS Comput. Biol. 11 e1004539
[56] Chen M M Guo D Q Wang T B Jing W Xia Y Xu P Luo C Valdessosa P A Yao D 2014 PLoS Comput. Biol. 10 e1003495
[57] Chen M M Guo D Q Yang X Yao D Z 2017 Frontiers in Computational Neuroscience 11 31
[58] Xia X F Wang J S 2014 Acta Phys. Sin. 63 140503
[59] Lytton W W 2008 Nat. Rev. Neurosci. 9 626
[60] Mccarthy M M Moorekochlacs C Gu X Boyden E S Han X Kopell N 2011 Proc. Natl. Acad. Sci. USA 108 11620
[61] Wang M L Wang J S 2015 Acta Phys. Sin. 64 108701
[62] Grant G Ernst N 2016 PLoS Comput. Biol. 12 e1005121